# Loading In Data and Libraries -------------------------------------------
setwd("")
data<-read.csv("02_Clean Racial Shift Replication.csv")

library(lavaan)
library(semTools)

# Fit Measures ------------------------------------------------------------
FIT <- c("chisq", "df", "cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")


# Treatment Check ---------------------------------------------------------
m_1d <- '
sc =~ SC_1 + SC_2 + SC_3 + ST_1 + ST_2 + ST_3 + PT_1 + PT_2 + PT_3 + RT_1 + RT_2 + RT_3
'
out_1d_geog <- cfa(m_1d, subset(data, Treatment == "Control"))
summary(out_1d_geog, fit.measures = T, standardized = T)

out_1d_race <- cfa(m_1d, subset(data, Treatment == "Population Shift"))
summary(out_1d_race, fit.measures = T, standardized = T)

out_1d_pol <- cfa(m_1d, subset(data, Treatment == "Political Shift"))
summary(out_1d_pol, fit.measures = T, standardized = T)

out_1d_cult <- cfa(m_1d, subset(data, Treatment == "Cultural Shift"))
summary(out_1d_cult, fit.measures = T, standardized = T)

m_4d <- '
sc =~ SC_1 + SC_2 + SC_3
st =~ ST_1 + ST_2 + ST_3
proto =~ PT_1 + PT_2 + PT_3
rt =~ RT_1 + RT_2 + RT_3
'
out_4d_geog <- cfa(m_4d, subset(data, Treatment == "Control"))
summary(out_4d_geog, fit.measures = T, standardized = T)

out_4d_race <- cfa(m_4d, subset(data, Treatment == "Population Shift"))
summary(out_4d_race, fit.measures = T, standardized = T)

out_4d_pol <- cfa(m_4d, subset(data, Treatment == "Political Shift"))
summary(out_4d_pol, fit.measures = T, standardized = T)

out_4d_cult <- cfa(m_4d, subset(data, Treatment == "Cultural Shift"))
summary(out_4d_cult, fit.measures = T, standardized = T)

# Output for table
summary(compareFit(out_1d_geog, out_4d_geog))
summary(compareFit(out_1d_race, out_4d_race))
summary(compareFit(out_1d_pol, out_4d_pol))
summary(compareFit(out_1d_cult, out_4d_cult))


# * ME --------------------------------------------------------------------
fit_config <- cfa(m_4d, data, group = "Treatment", 
                  group.equal = "none")
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[FIT]
# modindices(fit_config)

fit_metric <- cfa(m_4d, data, group = "Treatment", 
                  group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[FIT]
# modindices(fit_metric)

fit_scalar <- cfa(m_4d, data, group = "Treatment",  
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[FIT]
# modindices(fit_scalar)

comb_fit_all <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit_all

pchisq(comb_fit_all[3, 1]-comb_fit_all[1, 1],
       comb_fit_all[3, 2]-comb_fit_all[1, 2],
       lower.tail = F)

# output for table
summary(compareFit(fit_config, fit_scalar))
